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12 Abstract 

13 A new model coupling scheme with remote sensing data assimilation was developed for 

14 estimation of daily actual evapotranspiration (ET). The scheme represents a mix of the VegET, a 

15 physically based model to estimate ET from a water balance, and an event driven phenology 

16 model (EDPM), where the EDPM is an empirically derived crop specific model capable of 

17 producing seasonal trajectories of canopy attributes. In this experiment, the scheme was 

18 deployed in a spatially explicit manner within the croplands of the Northern Great Plains. The 

19 evaluation was carried out using 2007-2009 land surface forcing data from the North American 

20 Land Data Assimilation System (NLDAS) and crop maps derived from remotely sensed data of 

21 NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS). We compared the canopy 
parameters produced by the phenology model with normalized difference vegetation index 
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23 (NDVI) data derived from the MODIS nadir bi-directional reflectance distribution function 

24 (BRDF) adjusted reflectance (NBAR) product. The expectations of the EDPM perfonnance in 

25 prognostic mode were met, producing detennination coefficient (r ) of 0.8 ±0.15. Model 

26 estimates of NDVI yielded root mean square error (RMSE) of 0.1 ±0.035 for the entire study 

27 area. Retrospective correction of canopy dynamics with MODIS NDVI brought the errors down 

28 to just below 10% of observed data range. The ET estimates produced by the coupled scheme 

29 were compared with ones from the MODIS land product suite. The expected r =0.7 ±0.15 and 

30 RMSE = 11.2 ±4 mm per 8 days were met and even exceeded by the coupling scheme 

31 functioning in both prognostic and retrospective modes. Minor setbacks of the EDPM and 

32 VegET performance (r about 0.5 and additional 30 % of RMSR) were found on the peripheries 

33 of the study area and attributed to the insufficient EDPM training and to spatially varying 

34 accuracy of crop maps. Overall the experiment provided sufficient evidence of soundness and 

35 robustness of the EDPM and VegET coupling scheme, assuring its potential for spatially explicit 

36 applications. 

37 1. Introduction. 

38 There is growing consensus in the climate science community that the ability to precisely 

39 partition energy and matter fluxes on the land surface is key to improving our understanding of 

40 mesoscale atmospheric dynamics, ecosystem, responses to climate change, and interactions with 

41 human activities and institutions [ Pitman , 2003; Ibanez et al., 2010]. Since Manabe [1969], land 

42 surface modules (LSM) have become increasingly complex modules within most general 

43 circulation models (GCMs). The complexities of LSMs have grown substantially as scientists 

44 ask questions about the pace and consequences of climate change that require more precise 
answers. In pursuit of these answers, researchers have been coupling global and regional climate 
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models with a spectrum of modules detailing interactions between the land surface and the 

47 lowest level of the atmospheric boundary layer. Modules range from a set of simplified surface 

48 energy and water balance procedures to more detailed interactive systems like dynamic soil and 

49 vegetation modules, complete with light transfer, photosynthesis, and hydrological schemes. 

50 Computational resources often limit the level of detail in LSMs especially in regional studies that 

51 require finer spatial resolution. Also, there is a trade-off between the number of land surface 

52 characteristics that can be tracked and the greater spatial detail often needed for regional to local 

53 projects [Stensrud, 2007]. 

54 Applications of land surface models in regional studies were often focused on just a few 

55 variables of interest. In many instances, this narrower focus has led to the use of simplified 

56 schemes of land surface processes. Numerous local impact studies are turning to empirical 

57 methods based on relationships of modeled land surface characteristics to net radiation, 

58 precipitation, air temperature and other variables [. Nagler et al., 2005; Godfrey et al., 2007; 

59 Senay et al., 2007; Abramowitz et al., 2008; Jang et al., 2009; Gao et al., 2010], However, being 

60 developed on microclimatological data, empirical models were often unable to predict well when 

61 transferred to a different location, even under similar conditions \Li et al., 2009]. Yet, 

62 deployment of process-based LSMs to address local questions are often hindered by 

63 computational expense and a lack of appropriate ground level data to calibrate and validate at 

64 the level of spatial detail required. Also, several studies have expressed concerns about model 

65 assumptions, process parameterizations, and a limited range of parameters available for tuning 

66 [Sabater et al., 2007; Kiniry et al., 2008; Kang et al., 2009; Stancalie et al., 2010], all of which 

67 increase doubts about the likelihood of successful deployment of LSMs in regional to local 
studies. Alternatives solutions are needed to provide robust schemes capable of replacing 


68 



69 complex LSMs in finer spatial resolution studies. This paper presents a recent development in 

70 land surface modeling combining both physics-based and empirical approaches to take 

71 advantage of the strengths of each approach while yielding results on an appropriately fine scale. 

72 Our research focuses on how potential futures for rainfed agricultural production in the Northern 

73 Great Plains may affect regional hydrometeorology. Actual evapotranspiration (ET a ) was the key 

74 flux of interest. We chose to use a simplified simulator of ET a called VegET [Senay 2008], 

75 Similar to Godfrey et al. [2007], Kang et al. [2009] and Yuan et al. [2010], Senay ’s scheme relies 

76 on the Penman-Monteith equation [Monteith, 1964] to calculate reference ET (ETo) and handles 

77 the influences of soil water status and canopy phenology through the two coefficients: K s for 

78 soil water status and Kc P , for canopy phenology. The Penman-Monteith method is a physics- 

79 based one source model of evapotranspiration in cereal crops with fully developed canopies, 

80 used extensively by FAO [ Allen et al., 1998]. A key innovation of VegET is the modulation of 

81 ET a by a canopy phenology coefficient using a climatology of the normalized difference 

82 vegetation index (NDVI) as observed from spacebome sensors [Senay, 2008], 

83 The original implementation of VegET, however, could not serve our purpose because we were 

84 seeking how ET a would change in response to both interannual variability and changes in the 

85 crop area. Since they were derived from averages of past observations, a static retrospective 

86 climatology for K cp would not reflect changes in growing conditions [ Godfrey et al., 2007; 

87 Wegehenkel, 2009] or in the extent of cultivation [Kovalskyy and Henebry, 2011b]. Therefore, 

88 we replaced a static phenological parameterization with an interactive vegetation growth module. 

89 The use of fully functional crop growing modules with energy balance models in point based 

90 studies has been common practice [Maruyama and Kuwagata, 2010; Sancalie et al., 2010]. 
However, our study case required spatially explicit ET a estimates that would entail additional 
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92 parameterization, tuning and running time for the models like ALMANAC \Kiniry et al., 2008], 

93 CERES [ Mearns et al., 1999], CROPWAT [ Sancalie et al., 2010] or MODWht [Kang et al., 

94 2009]. Moreover, these specific crop models did not have freely available versions capable of 

95 working with raster inputs and producing spatially explicit estimates. Conversely, the vegetation 

96 growth modules in global LSMs were developed to deliver spatially explicit results [Dickinson et 

97 al., 1998; Foley et al., 2000; Bondeau et al., 2007; Campo et al., 2009]. Even the most advanced 

98 modules do not provide crop specific canopy behavior; instead, they were designed to mimic 

99 seasonal patterns of very broad classes of vegetation functional types [ Bonan et al., 2003; 

100 Lawrence and Chase, 2007]. 

101 Here we have used the Event Driven Phenology Model (EDPM), which was recently developed 

102 as a phenology model that can simulate seasonal dynamics of canopy properties (e.g., in terms of 

103 a vegetation index) [Kovalskyy and Henebry, 2011a, 2011b]. The model was shown to capture 

104 fine temporal details of canopy behavior [Kovalskyy and Henebry, 2011a, 2011b] which has 

105 been called “crucial” for ET and other surface fluxes [ Dickinson et al. 1998; Foley et al. 2000; 

106 Pitman, 2003; Gorfrey et al. 2007; Prihodko et al. 2008; Rosero et al. 2009; Rotzer et al. 2010; 

107 Zha et al. 2010], The EDPM uses virtually the same set of forcings as the Penman-Monteith 

108 equation to build seasonal trajectories of canopy properties. Essentially, the model provides a 

109 computationally inexpensive replacement for a dynamic vegetation model with a phenology sub- 

110 module. The model also has an option of a simple, fast ID data assimilation scheme for satellite 

111 observations which is a great advantage for spatially explicit simulation studies. The EDPM has 

112 been coupled with VegET and evaluated against flux tower observations of ET a [Kovalskyy and 

113 Henebry, 2011b], where it performed better or comparable to the results obtained by Nagler et 

114 al. [2005] and Abramowitz et al. [2008], 



115 This paper presents an assessment of the perfonnance of the EDPM on its own and also in 

116 conjunction with VegET within a spatially explicit application. Our task was to select 

117 appropriate sources of scientifically sound data products that would enable pixelwise 

118 comparisons of daily canopy states and ET a estimates. We were looking to assess two aspects of 

119 the coupled model perfonnance. First and foremost, we focused on temporal and spatial 

120 behavior of differences between our estimates and reference data. Analyzing the results from the 

121 three years (2007-2009) within the study region delimited by croplands of Northern Great Plains, 

122 we tried to capture both inter-annual and intra-annual variability of residuals as well as 

123 conelation between reference data and estimates produced by our model. Second, we looked at 

124 the ability of the EDPM to capture the key dates of the three growing seasons. We contrasted 

125 phenological dates reported to National Agricultural Statistics Service (NASS) by fanners with 

126 the dates produced by the EDPMs phenophase control module. Specifically, the Start of Season 

127 (SoS), End of Season (EoS), and Length of Season (LoS) became the main criteria for the 

128 evaluation. We also tried to incorporate spatial and temporal variability of phenological metrics 

129 into the evaluation process. This assessment helped us to identify the strong points of the EDPM 

130 and to prioritize directions for model improvement. 

131 2. Methods and materials, 

132 2.1 Study area. 

133 The study area includes Nebraska, Iowa, Minnesota, North and South Dakota entirely and parts 

134 of Illinois and Indiana. Together these states have more than half of the nation’s maize and 

135 soybean crops and comprise the major part of the US corn and soybean belts. There strong 

136 gradients of ET across the region. The northern tier has only 600 mm ET annually; whereas at 

137 the southern end, the annual ET can reach 1000 mm, but only 400 mm at the western extreme. 
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Maize and soybean are the most prevalent crops across the region. Farmers use different genetic 
varieties of these crops to match the growing conditions of their farms [Ransom et al., 2004]. 
The green-up of the area starts in early May on the southeast end but for the northwestern part of 
the region it can happen as late as mid- June if spring comes late. The length of the growing cycle 
also varies greatly; it can last almost five months in the south and barely more than three months 
in the north. 



Figure 1. The study area (dark gray) depicted as at least 50% corn/soybean crop cover 
during 2007-2009. 

2.2 Coupling the VegET and the EDPM in a spatially explicit manner. 

The idea motivating the development of VegET was the use the time series of remotely sensed 
vegetation indices to drive the canopy factor that modifies ETo as calculated by the Penman- 
Monteith model. The original design of VegET [Senay, 2008] used very simple empirical 
transformations from the normalized difference vegetation index (NDVI) to phenology driven 
coefficients based on thresholds and the observed variability in NDVI climatologies derived 
from long-tenn AVHRR observations. However, Kovalskyy and Henebry [2011b] demonstrated 



154 the use of an interactive event driven phenology model [ Kovalskyy and Henebry, 2011a] to 

155 replace the climatologies in the VegET for point-based estimation of daily ET a . The coupling 

156 scheme was shown to account for contemporaneous fluctuations in the canopy component of 

157 evapotranspiration 

158 The experiment described here evaluates the performance of the EDPM and the VegET after the 

159 coupled models were extended for deployment in a spatially explicit manner. In addition to 

160 simulating the temporal dynamics of maize and soybean over the three year period (2007-2009), 

161 the EDPM was also providing seasonal canopy trajectories for a third vegetation type: grassland. 

162 Using weather forcing from the North American Land Data Assimilation System (NLDAS), the 

163 EDPM transfonned the data in to events (rain, heat stress, frost, insufficient insolation, adequate 

164 insolation, and growing degree-days) and further produced daily values of Tower NDVI 

165 [Huemmrich et al, 1999], The model simultaneously estimated phenological transition dates in 

166 the three growing seasons for each vegetation type (Fig. 2). The TNDVI trajectories were mixed 

167 linearly based on the proportion of their cover within areal units (0.05 degree pixels) and later 

168 transformed into phenology coefficients as described in Kovalskyy and Henebry [2011b]. The 

169 percentages of cover for each crop and grassland were taken from MODIS based crop maps 

170 products [ Chen et al., 2007] and aggregated into standard 0.05 degree (-5 km) pixels that form 
the spatial unit of analysis for this investigation. 
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Figure 2. Data processing scheme for the experiment. Rounded boxes are modeling and 
data preparation procedures; stacks are image time series (ITS) of data; squared boxes are 
maps of results. 

Using the workflow shown in Figure 2 the coupled scheme was tested in the prognostic mode 
(running the forcing only) and the diagnostic mode (involving data assimilation scheme with 
MODIS NDVI observations). The use of data assimilation techniques is becoming increasingly 
popular in evapotranspiration studies [ Meng et al., 2009; Anderson et al., 2010; Miralles et al., 
2010; Godfrey and Stensrud, 2010]. Most of these projects relied on remotely sensed data to 
improve their estimates of ET addressing the spatial variability of land surface. While bringing 
improvements in performance, these techniques have been criticized as being temporally 
constrained and scene dependent [Li et al., 2010]. Our study took a more general approach to 
data assimilation using an unambiguous relationship between Tower NDVI and MODIS NDVI 


186 established previously \Kovalskyy et al., 2011]. Relying on this relationship Kovalskyy and 

187 Henebry [2011a] presented a one-dimensional Kalman filter (1DKF) data assimilation scheme in 

188 which the EDPM used MODIS NBAR NDVI to adjust its estimates of canopy states. 

189 2.2. 1 Differences from point based deployment of the EDPM. 

190 Several features in the EDPM model were added and modified so that the model could represent 

191 spatial variability of canopy development during growing season. First of ah, the model received 

192 the ability to represent pixels with mixed vegetation cover. The Figure 2 shows the linear mixing 

193 procedure used to derive values of vegetation index in a pixel with partial covers of grassland 

194 and the two crops. Direct linear mixing of NDVI values based on their share of pixel area has 

195 been criticized in the literature due to its impact on outcomes [ Settle and Campbell, 1998; Roy, 

196 2000; Busetto et al., 2008]. Since the NDVI is not a linear function of red and near infrared 

197 reflectances, linear mixing should be performed on reflectances first so that “unbiased” NDVI 

198 values can be obtained later. However, a relatively small impact from direct linear mixing 

199 (DLM) of NDVI values may not be entirely prohibitive since the reflectances coming from 

200 MODIS products have their own errors of the estimate [ Roy et al., 2005], If the differences 

201 between the DLM NDVI and the NDVI derived from reflectances can be kept within the margin 

202 of error propagated into the unbiased NDVI, then one can successfully perfonn linear mixing 

203 directly using NDVI values. The magnitude of differences with true values depends on the 

204 number of endmembers used in linear mixing and varies greatly across space and time. We 

205 evaluated a real (empirical) example to demonstrate how the direct linear mixing procedure 

206 impacts the resulting values of NDVI. 

207 First we took a 1000 by 1000 pixel subset from the MODIS Nadir Bi-directional Reflectance 

208 Distribution Function (BRDF) Adjusted Reflectance (NBAR) MOD43A4 [Schaaf et al., 2002] 



209 version product covering the central part of the study area. We screened for snow and clouds and 

210 based on averages aggregated the 500m reflectance values from MODIS bands 1 and 2 to 

211 produce image time series with 1000 m by 1000 m, 2500 m by 2500 m, 2500 m by 5000 m 

212 (rectangular shape pixel) and 5000 m by 5000 m pixel sizes. For each time series less than 25 

213 km“ in size, we calculated NDVI [1] that was later mixed directly into 5000 m pixels. 

214 Correspondingly, the four sets of results represented 100, 25, 4, and 2 endmember mixing. Out of 

215 5000 m reflectance data we calculated “true” NDVI [1] and expected error [2] propagated from 

216 reflectances: 0.004 for band 1 and 0.015 for band 2 [ Roy et al., 2005]. 
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219 where p N and p R are the reflectance values of near infrared and red bands respectively, and a“ N 

220 and a r are the associated variances. In this setup where all resolutions of NDVI data were nested 

221 within 5000 m pixels, we expected to see the difference coming just from linear mixing without 

222 other effects such as re-projection or resampling that may otherwise contribute to the difference 

223 [Roy, 2000], Figure 3 demonstrates the temporal dynamics of average differences between true 

224 NDVI and the four DLM NDVI sets coming from linear mixing with different number of 

225 endmembers. 
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Figure 3. Consequences of linear mixing - an observation based example. Difference NDVI = Mixed 
NDVI minus “unbiased” NDVI. 

It is seen clearly from the figure above that the impacts (differences) from linear mixing increase 
as the number of endmembers grows. The differences become significant when the number of 
mixing endmembers reaches 4. In our study we used only 3 endmembers (maize, soybeans and 
grassland) to be mixed into 0.05 by 0.05 degree pixel representing the trajectory of Tower NDVI 
values. Considering the magnitude of errors from demonstrated direct linear mixing examples it 
was safe for us to assume that 3 endmember linear mixing did not make a significant impact on 
the seasonal trajectories of TNDVI produced by the EDPM. We also have to point out here that 
the minor impacts from linear mixing appeared to be negligible (20 times smaller) compare to 
the estimate errors of the EDPM reported in Kovalskyv and Henebry [2011a]. In this context, the 
impacts from direct linear mixing could hardly make a difference for comparisons undertaken in 
this experiment. 


240 Next, the transformation of mixed TNDVI into phenology driven coefficient K cp had to be 

241 generalized (unified). In the prior point based studies we found TNDVI and Kc P relationships to 

242 be different for crops and grassland \Kovalskyy and Henebry 2011b]. For the later land cover 

243 type, the linear model carried substantial noise that we tried to compensate with modeling of 

244 residuals through their relationship with vapor pressure deficit. We did not find the same 

245 relationship in residuals for crops assuming the bias in grassland was due to differences in 

246 equipment calibration. Therefore, we used a single linear model with the slope of 1 .22 and offset 

247 of 0.01 to transform modeled TNDVI into K cp in this spatially explicit experiment. This 

248 relationship was derived on observations of TNDVI and K cp on crops and proved its efficacy in 

249 Kovalskyy and Henebry [2011b]. 

250 Finally, the spatially explicit application of the event driven phenology model required some 

251 amendments in the functioning of the phenological phase control module described in Kovalskyy 

252 and Henebry [2011a]. Trained on specific locations and tested on locations with similar climatic 

253 conditions, the EDPM required a supplementary mechanism to match the variability of the 

254 growing season dates within a much wider range of conditions than in the initial testing studies. 

255 Latitudinal gradients for the emergence of vegetated cover which marks the start of the season 

256 were applied to the two controlling variables: thennal time and elapsed days since January 1. 

257 For the elapsed days we applied a 4 days per degree northward gradient suggested by Hopkins 

258 [1918]. The thermal time triggers for the onset of greening was also modified with the latitude 

259 using slopes and intercepts tuned for three vegetation types following the phenological transfer 

260 functions described in Henebry [2010], To deal with the southward increase in the duration of 

261 the growing seasons, we adjusted the dynamic triggering for transitions between phenological 

262 phases. The adjustment made the transition probability threshold vary inversely with the latitude. 



263 This helped to postpone transitions between phenological phases for locations to the South of 

264 training sites, while accelerating the transitions to the North. 

265 2.3 Data sources and preparations for the experiment 

266 The experiments conducted within this investigation had to use various data sources to reach 

267 their goals: (1) running the EDPM plus VegET coupling scheme required weather forcing data; 

268 (2) percent crop cover data were necessary for the EDPM to produce seasonal canopy trajectories 

269 of pixels with mixed vegetation cover; (3) NDVI observations were needed to verify the 

270 EDPM’s prognosis of seasonal canopy trajectories and later to produce retrospective outcomes; 

271 (4) observations of actual ET were needed to evaluate the quality of estimates produced by the 

272 EDPM plus VegET coupling scheme; and (5) crop progress reports were crucial for assessment 

273 of the EDPM module responsible for estimating dates of phenological transitions. 

274 The meteorological forcings for the EDPM and the VegET were supplied by the North American 

275 Land Data Assimilation System (NLDAS) in native GRIB1 format (1 hour temporal and 0.125 

276 degree spatial resolutions). The choice of NLDAS [ Mitchell et al., 2004] was based in part on the 

277 fact that these forcings were validated on the southern Great Plains adjacent to our study region 

278 \Luo et al., 2003]. The original time series of weather data were aggregated into daily image 

279 time series and resampled into 0.05 degree grid using nearest neighbor procedure to match the 

280 MODIS Climate Modeling Grid (CMG) projection. The last transformation preserved most of 

281 the original data and allowed for the fusion of MODIS NDVI data with the EDPM produced 

282 seasonal canopy trajectories and later calculation of the ET a at 0.05 degree resolution. The list of 

283 forcing variables included: 2 meter air temperature [K]; 2 meter specific humidity [kg/kg]; 

284 surface pressure [Pa]; U wind component [m/s]; V wind component [m/s]; downward shortwave 

285 radiation [W/m“]; downward longwave radiation [W/nT]; total precipitation [kg/m ]. The forcing 



286 dataset and LSM simulations of NASA’s Mosaic model from NLDAS Phase 2 were obtained 

287 from the NASA Goddard Earth Sciences Data and Information Services Center (GES DISC) at: 

288 http://disc.sci.gsfc.nasa.gov/hydroIogy/data-hoIdings. 

289 Land cover data came in the fonn of MODIS based crop maps for 2007 through 2009. The 0.5 

290 km resolution maps were provided directly by members of the product development team 

291 [ Chang et al, 2007]. The procedures for deriving the percent covers of maize and soybeans were 

292 based on decision tree techniques applied on level 2 MODIS reflectance in seven bands covering 

293 visible and infrared portions of the electromagnetic spectrum. Additional metrics capturing the 

294 temporal development of land surface properties relevant to the vegetation were also fused into 

295 the procedure. The 2007 paper reported drawbacks in using the universal sampling approach in 

296 the development of their decision tree model which were related mostly to the differences in 

297 cropping timing and radiometric properties of underlying soils in different areas of CONUS. As 

298 an alternative, they proposed a single state based modeling of percent crop type cover which 

299 significantly improved the performance of the procedure in independent tests. 

300 Here we were able to use the latest versions of crop cover maps derived from state based 

301 decision tree models of Nebraska, Iowa, Minnesota, North and South Dakota, where the most 

302 variability was captured in the least amount of training. Some adjacent areas of other states were 

303 combined in the final areal compositing procedure. Within delineated study area where each 

304 pixel had at least 50% crop cover, the proportion of grassland was assumed to be the remainder 

305 of a pixel cover. This assumption was based on the NLDAS land cover scheme that considered 

306 grassland as the second most abundant land cover within our region [Luo et al. , 2003], We also 

307 masked out all non-grassland or non-cropland land covers from our study area based on the 



308 MODIS land cover product (MCD12C1, IGBP classification type available at 

309 ftp://e4ftl0 1 .cr.usgs.gov/MOTA/MCD 12C 1 .005). 

310 Verifying the estimated canopy states and actual evapotranspiration on a spatially explicit basis 

311 posed some difficulties only because there are so few observational datasets available for such 

312 analyses. One such comparison was the TNDVI time series generated by the EDPM in 

313 prognostic mode with the MODIS NDVI time series. The reference NDVI image time series 

314 were produced from MODIS NBAR data (MCD43C4 version) available at 

315 ftp://e4ftl01.cr.usgs.gov/MOTA/MCD43C4.005. First, bands 1 and 2 of NBAR data in CMG 

316 projection were extracted and screened for insufficient quality records using QA bits. Then, the 

317 NDVI [1] was computed out of screened red and near infrared reflectance data and organized 

318 into image time series. 

319 Potentially, daily ET a estimates from the EDPM plus VegET scheme had several sources of 

320 reference data since at the time of our study two ET monitoring products were on their way to 

321 public release [Mu et al, 2007; Anderson et al, 2010]. However, only the MODIS 

322 evapotranspiration product (MOD16) data had become publically accessible [Mu et al., 2009]. 

323 Therefore the MOD 16 product had become our first choice for reference when assessing the 

324 quality of results from the coupled EDPM+VegET scheme. This product presents estimates of 8 

325 day sums of actual and reference ET modeled from weather forcings and remotely sensed 

326 properties of the land surface [Mu et al., 2007]. Standard HDF files were obtained from 

327 ftp.ntsg.umt.edu/pub/MODIS/Mirror/MOD16/MOD16A2.105_MERRAGMAO/. The original 

328 actual ET layers of MOD16A2 version of the product with 1 km resolution were spatially 

329 aggregated and then resampled into 0.05 degree grid again to match the MODIS CMG projection 

330 adopted as the basis for this experiment. 



331 The M0D16 product has been closely approaching the in situ measured ET a with each 

332 improvement to its procedures [Mu et al., 2009, 2011], yet it is still a product with varying 

333 degree of spatial and temporal uncertainty. Therefore, we retained an alternative set of ET a 

334 estimates with which to compare our results. We selected the outcomes of NASA’s Mosaic LSM 

335 [ Koster and Suarez, 1994, 1996] from NLDAS as an alternative reference point for comparison 

336 based on the validation studies of Mosaic LSM [Koster and Suarez, 2003; Koster et al., 2004], 

337 To match the fonnatting of the first reference product (MOD 16), the daily ET a estimates from the 

338 coupled EDPM + VegET scheme and from Mosaic LSM were each temporally aggregated into 8 

339 day ET a totals. 

340 The accuracy in estimating phenological dates has always been a subject of point based 

341 validation studies [Menzel et al., 2006; Schwartz et al., 2006; Richardson et al., 2009; Zhang et 

342 al., 2009; White et al., 2009; Dufour and Morin, 2010]. In the search of reference data we 

343 examined the National Agricultural Statistics Service (NASS) weekly Crop Progress reports on 

344 the percentages of crops achieving crop specific phenophases. From this source we could only 

345 obtain information about growing season progress on the state level since the county level 

346 reports were inconsistent. Therefore, we reorganized the pixel based EDPM reports into the daily 

347 state level growing season progression time series to see the parallels between reported and 

348 observed dates. These dates were compared with two available years (2008, 2009) of state level 

349 crop progress reports obtained for the five states (Nebraska, Iowa, Minnesota, North and South 

350 Dakota) from the NASS archives: 

351 http://www.nass.usda.gOv/Data_and_Statistics/Quick_Stats_l.0/index.asp 

352 Considering the spatial mismatch, temporal precision differences, and the differences in 

353 biogeophysical meaning between reported events and dates estimated by the EDPM, we have 



354 chosen to rely mostly on the midpoints of distributions in phenological metrics for our 

355 comparison. Therefore in the analysis we used midpoint dates (when 50% of crops went through 

356 start of season [SoS] or end of season [EoS] ) and their inter-quartile range (IRQ) as a measure of 

357 data variability. Based on SoS and EoS dates we also calculated the lengths of seasons (LOS) 

358 together with their inter-quartile ranges. The LoS values from the NASS reports were calculated 

359 by subtracting the 50% EoS date from 50% EoS date and the IQR 75% EoS date minus 25% EoS 

360 date and 25% EoS date minus 75% EoS date. The IQR in the LoS data generated during our 

361 experiment were collected directly from the EDPM reported pixel phenology dates. 

362 2.4 Road map for analysis. 

363 Resulting test runs of the EDPM and the coupled scheme with VegET produced several sets of 

364 results for the evaluation. First, the image time series of TNDVI estimated by the EDPM in 

365 prognostic mode were compared with MODIS NBAR NDVI data. Despite the discrepancy in 

366 temporal resolution (8 day for MODIS products and daily for our estimates), the comparison 

367 could give a good idea of how close our predictions were to the observations. In preparation for 

368 such comparison, the EDPM outcomes went through the transfonnation into MODIS NDVI 

369 using the relationship developed in Kovalskyy and Henebry [2011a] and confirmed in Kovalskyy 

370 et al. [2011]. Avoiding the comparison of data beyond the growing season where the EDPM 

371 cannot produce TNDVI, we allocated only the results and reference data representing the period 

372 from early March (97 th day of the year) to late October (305 th day of the year). In addition to that, 

373 only the dates matching the beginnings of 8 day compositing periods of MODIS products (not 

374 the averages over compositing period) were selected for comparison. 

375 In diagnostic mode the EDPM used the former reference— MODIS NBAR NDVI data — to 

376 correct its outcomes via the built-in data assimilation scheme [ Kovalskyy and Henebry , 2011a]. 



377 Therefore, to assess the model performance in diagnostic mode, we had to rely on error 

378 propagation to infer the accuracy of the assimilation-enhanced EDPM estimates of TNDVI. 

379 Prognostic and diagnostic versions of the EDPM outcomes were used to parameterize VegET to 

380 produce corresponding ET a outcomes. Aggregated into 8 day totals to match the fonnat of first 

381 reference data, the ET a estimates from both prognostic and diagnostic runs of the scheme were 

382 compared with the temporally matching image time series of actual evapotranspiration from 

383 MOD 16 product validated by Mu et al. [2009] and Mosaic LSM validated by Luo et al. [2003]. 

384 Only the time series of ET a from early March to late October were used in the comparisons. 

385 In our assessment we relied generally on the two most common measures of perfonnance: 

386 coefficient of determination (r~) and root mean square error (RMSE). The first measure showed 

387 the ability of produced estimates to follow the observed developments of the modeled variable. 

388 RMSE showed the overall level of departure of modeled TNDVI and ET a from what we assumed 

389 to be the reality (reference datasets). Based on the results received in Kovalskyy and Henebry 

390 [2011a] and Kovalskyy and Henebry [2011b], the expected perfonnance levels for the canopy 

391 state estimates (viz., NDVI) simulated by the EDPM were r 2 =0.8 ±0.1 and RMSE = 0.1 ±0.025. 

392 For ET a , the expected performance levels were r =0.7 ±0.15 and RMSE =1.4 ±0.5 mm per day, 

393 but transformed into 8 day values by simple multiplication yields RMSE = 11.2 ±4 mm per 8 

394 days. Additionally, the results were examined for the presence of biases in the residuals. 

395 Analyzing differences with reference data, we aimed to assess both temporal and spatial aspects 

396 of their distributions to receive clear contrasts between sets of our modeling results and reference 

397 data. 

398 In its collection, the NASS archive offered emergence and maturity dates for maize as well as 

399 emergence and leaf drop dates for soybeans. We assumed these phenological turn points to be 
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closely related to the SoS and EoS dates produced by the EDPM. Comparing phenological data 
we plotted our estimates against references expecting to see connections between plant 
physiological events and their manifestation in the temporal dynamics of optical properties of the 
vegetated surface. 

3. Results. 

3.1 Contrasting the EDPM derived NDVI against MODIS product. 

The maps representing performance measures for each year were produced to show how the 
ability of the EDPM to represent the canopy conditions varies in space. We also included the 
maps of average seasonal propagated errors into Figure 4 from results received after the data 
assimilation (retrospective mode) to contrast those with RMSE obtained during uncorrected 
(prognostic) estimation. 


(a) (b) (c) 



0 . 0 - 0.5 0 . 5 - 0.8 0 . 8 - 1.0 0.12 0.24 0.36 0.48 0.12 0.24 0.36 0.48 



412 


Figure 4. Comparison of the EDPM produced vegetation index against MODIS NDVI 

413 within the study area, («) Coefficient of determination (r 2 ); ( b ) Root mean square error; (c) 

414 Seasonally averaged propagated daily NDVI error after assimilation of MODIS NDVI 

415 observations. 

416 The figure above clearly demonstrates that the EDPM was well fit for to the task of following the 

417 dynamics of observed MODIS NDVI. Maps in the left column are dominated by dark color 

418 representing r of 0.8 and higher. The r values had a tendency to decrease toward the borders of 

419 the study area and whereas 2007 was the year with the worst performance, 2008 the best. The 

420 same conclusion was supported by the RMSE maps in Figure 4. The overall level of error 

421 reached 0.18 for 2007, but dropped to just above 0.11 for 2008. The right column of Figure 4 

422 shows the unifonn distribution of average seasonal propagated errors throughout the study area 

423 after EDPM predictions were updated with MODIS NDVI observations. The general level of 
propagated errors was very close for all three years and constituted slightly less than 0.1. 
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Figure 5. Spatial distributions of residuals (NDVIedpm - NDVImodis): («) seasonal means of 
residuals; ( b ) standard deviations of residuals. 

Figure 5 above reveals that the EDPM was mostly underestimating the value of NDVI. Again the 
picture changed for different years and the character of bias reversed towards the peripheral areas 
of the study region. The year of 2007 came out as the most biased having the mean of residuals - 
0.2 to -0.3 spread along the western Iowa and Minnesota borders. For 2008 and 2009, most of 
the seasonally averaged differences between observed and modeled NDVI varied between -0.2 
and 0.1. The variability of the residuals grew from the center towards the borders for each year. 
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However, similar absolute values of RMSE (Fig. 4b) and mean residuals (Fig. 5a) point that the 
bias was rather uniform in time for most of the study area. 

A closer look into intra-annual dynamics of residuals (Fig.6) reveals similarities in developments 
seen in both the mean difference with observations and the standard deviation of residuals within 
the three growing seasons. 
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Figure 6. Temporal dynamics of residuals (NDVI edpm - NDVImodis) during the 2007-2009 
growing seasons. Light grey squares represent season of 2007; darker grey diamonds are 
2008; and black triangles are 2009. 

The trajectories in Figure 6 represent temporal dynamics of residuals averaged over the entire 
study region (18.7k pixels). It is seen clearly that the biases from different years went through 
similar seasonal patterns. The graphs show that the EDPM in prognostic mode was starting up 
seasons with minor underestimation and kept it at this level till the growing season started for 


447 maize and soybeans. The mean of residuals was dropping at every phenological transition point 

448 which constitutes a source of perfonnance problems in the EDPM [ Kovalskyy and Henebry, 

449 2011a, 2011b]. After the change of the phenological phase, the differences with observations 

450 came back to the initial level. This pattern means that corrections of the model outcomes during 

451 phase change were needed to decrease the bias and make the bias more stable. Overall, the 

452 analysis of the EDPM performance suggests that although the errors from EDPM were higher, 

453 they were still within the expected range based on prior performance. However, the results from 

454 the EDPM can be found reasonably accurate for prognosis or retrospective temporal gap filling 

455 in observations, considering the fact that the NDVI from EDPM carries the uncertainty from 

456 transformation to MODIS NDVI (standard error of the slope of 0.1 I [Kovalskyy and Henebry, 

457 2011a]), and the uncertainties of the crop maps proliferated through the mixing process. 

458 3.2 Contrasting ET„ estimates from the EDPM+VegET scheme against MODI 6 product. 

459 Before comparing the results from the EDPM+VegET with references, it is important to note that 

460 the gaps between the two reference datasets were substantial. Plot ( a ) in Figure 7 clearly shows 

461 that compared with MOD 16 product, Mosaic ET a first overestimated and then brought bias close 

462 to 0 in the middle of the growing season, but later it returned to overestimation. The two versions 

463 of the EDPM+VegET estimates representing ET a derived with and without assimilation via 

464 1DKF scheme also had their differences shown in plot ( b ) of the figure 7. Following the 

465 previously noted pattern of underestimation of canopy properties by the EDPM working in 

466 prognostic mode, the prognosis of ET a values was lower than ET a produced in diagnostic mode 

467 (with 1DKF). The variability of residuals in Figure 7b exhibited similar temporal behavior to the 

468 one found in the bottom plot of Figure 6. 
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Figure 7. Temporal dynamics of ET a residuals during the 2007-2009 growing seasons. ( a ) 

ET a Mosaic - ET a MOD16; ( b ) ET a EDPM +VegET~ ET a EDPM with lDKF+VegET ? (c) ET a EDPM +VegET - 
ET a Mosaic? iff) ET a EDPM +VegET “ ET a MOD16? (&) ET a EDPMwith 1DKF +VegET " ET a Mosaic? if) ET a 


473 edpm with lDKF+vegET - ET a modi6- Light grey squares represent season of 2007; darker grey 

474 diamonds are2008; and black triangles are 2009. 

475 Retaining the main features from plots a and b, the remaining graphics of Figure 7 show the 

476 temporal dynamics of differences between two reference datasets and the two sets of 8 day ET a 

477 estimates from the EDPM+VegET scheme. In prognostic mode the EDPM+VegET results were 

478 starting growing seasons with underestimation of 15 mm per 8 days compare to ET a produced by 

479 Mosaic. In the midseason the difference came close to zero, but later a smaller (~10 mm per 8 

480 days) underestimation prevailed again (Fig. 7c). Meanwhile compared to the ET a from MOD 16, 

481 the prognosis from the EDPM+VegET showed close to 0 difference for most of the season with 

482 slight overestimation in early June (up to 7 mm per 8 days) and underestimation of the same 

483 magnitude in late August (Fig. 7d). The variability of residuals for prognostic ET a estimates 

484 remained high and had a clear temporal pattern driven by phenology. 

485 The estimates of ET a obtained with the EDPM+VegET working in diagnostic mode (with 1DKF) 

486 exhibited similar behavior of residuals when compared to reference datasets. Differences with 

487 Mosaic were negative at the beginnings of growing seasons (Fig. 7e), but in the mid-season the 

488 curves drifted toward slight (up 5 mm per 8 days) overestimation which later changed back to the 

489 underestimation of 15 mm per 8 days again (Fig. 7e). Compared with MOD 16 the 

490 EDPM+VegET diagnostic estimates produced residuals that signal slight overestimation early in 

491 the growing season. Later, however, the residuals came close to 0 and remained there till the end 

492 of growing season indicating good match (Fig. 7f). The variability of residuals for retrospective/ 

493 diagnostic ET a estimates from EDPM+VegET dropped quite dramatically in both comparisons 

494 (Fig. 7e,f) showing the relative efficacy of data assimilation for this method of ET a estimation. 
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Overall, the EDPM+VegET scheme showed closer temporal resemblance with MOD 16 product 
and therefore further we present figures representing the spatial particularities of the coupled 
model performance compared to the MODIS product. (Analogous figures showing the 
comparison with Mosaic can be found in Appendix A.) 
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Figure 8. Comparison of MOD16 ET a with the ET a produced by EDPM+VegET working 
in (A) prognostic mode and (B) diagnostic mode involving 1DKF assimilation. (/) 
Coefficient of determination (r 2 ); (ii) Root mean square error (mm per 8 days). 


Both parts of figure 8 show that EDPM+VegET scheme was able to follow the dynamics of ET a 
in the reference dataset and produced high values of determination coefficient exceeding the 
expectations set in previous section. Average coefficient of determination was above 0.8 level 
for the scheme working in both prognostic and diagnostic modes. In 2008, however, the average 
value of r dropped to the expected 0.7 level for both versions of derived ET a (Fig. 8A and B). 
The distribution of r values within the study area was more even in the results from the coupled 


510 scheme working in diagnostic mode involving 1DKF assimilation with MODIS NDVI data (Fig. 

511 8 A). In both modes the EDPM+VegET scheme showed lower r in the western peripheral 

512 regions where the accuracy of crop cover maps was lower. Correspondingly, the RMSE values in 

513 those regions were higher especially in the results of the scheme working in prognostic mode. In 

514 the results from diagnostic mode RMSE had more unifonn distribution and constituted around 6 

515 mm per 8 days on average which is half of what was expected. The average RMSE for 

516 EDPM+VegET outcomes derived in prognostic mode was about 8 mm per 8 days. Transformed 

517 into corresponding units, this perfonnance would be comparable to Nagler et al. [2005] or 

518 Abramowitz et al. [2008], if the ET a data from MOD 16 product approximated the reality with the 

519 accuracy of flux tower instruments [Mu et al., 2009]. A point based flux tower validation study 

520 has shown that the scheme can approximate daily ET a in crops with similar accuracy [Kovalskyy 

521 and Henebry, 2011b]. 
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Figure 9. Spatial distributions of residuals (A) ET aED pM+ve g ET- ET a modi6 (B) ET aE DPMwith 
iDKF+vegET- ET a modi6- (0 annual mean of residuals (mm per 8 days) ; (**) standard 
deviation of residuals (mm per 8 days). 

The contrast between the two sets of ET a estimates from the EDPM+VegET scheme can be 
easily depicted from the Figure 9 (A and B). In the left column (z) of panel A of Figure 9, the 
prognoses of ET a had mostly negative bias changing to overestimation in the peripheral areas of 
the study region (both east and west). The magnitudes of the mean residuals deviated not too far 
from 0 giving a peak of up to 12 mm per 8 days in 2007 in the central part of the study region. 
Left column (zz) of Figure 9A shows uneven distribution of variability in residuals revealing 



533 clusters of instability in perfonnance coming from EDPM+VegET scheme working in prognostic 

534 mode. Panel B of the Figure 9 shows that performance of the EDPM+VegET scheme was more 

535 stable during the work in diagnostic mode. The bias in the left column of the Figure 9A was 

536 mostly positive fluctuating no more than 9 mm per 8 days. There was less contrast between years 

537 and also less difference between various parts of the study region. Smaller and more 

538 homogenously distributed standard deviations of residuals (Fig. 9B column ii) also indicated a 

539 greater stability in performance compare to prognostic mode (Fig. 9A column ii). 

540 Contrasted with the ET a estimates from Mosaic (Appendix A) the results from the 

541 EDPM+VegET scheme were less correlated and had greater spatial variability in RMSE and 

542 residuals. Figures in Appendix A clearly demonstrate the problem in the central part of the study 

543 area (especially during 2007 growing season) that came from numerous differences in 

544 approaches to the ET a modeling and the associated assumptions made about the parameter 

545 datasets e.g. land cover types, soil types, LAI, etc [Koster and Suarez, 1996; Mitchell et al., 

546 2004], Nevertheless, the expected performance of r 2 =0.7 ±0.15 and RMSE = 11.2+4 mm per 8 

547 days were achieved by the coupled models working only in retrospective/diagnostic mode using 

548 MODIS observations for correction of simulated TNDVI trajectories. 

549 3.3 Comparison of growing season parameters. 

550 The need to evaluate the performance of the phenological control module in the EDPM was well 

551 motivated by the patterns in residuals seen in Figures 6 and 7. Therefore, we highlight the 

552 contrasts between the EDPM estimated and in situ start of season [SoS] and end of season [EoS] 


553 


dates reported to NASS. 
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Figure 10. Contrasting start and end dates of the growing season for the two crops and two 
years. 

Figure 10 shows fairly good agreement between observed and estimated parameters of the two 
growing seasons. It also reveals the persisting delays in SoS for maize crops within all five 
states. Nevertheless, the 2 weeks delays in SoS prognoses were comparable with errors 
encountered in retrospective analyses by Fisher et al. [2006] Zhang et al. [2009] and Kovalskyy 
et al. [2011]. Meanwhile, the estimates of both SoS and EoS for soybeans were even more 
precise and consistent. Figure 10, however, does not show the variability of the start and end of 
season dates where dramatic differences arise between NASS reports and the EDPM. To 
condense the graphical information, we brought the variability measure, the interquartile range 
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(IQR) into Figure 11, which also shows the scatterplots in length of season (LoS). Similar 
patterns occur in the variability of SoS and EoS (data not shown). 
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Figure 11. Contrasts between the EDPM estimates and NASS reports in the length of 
season and its variability for the two crops and two years. 

The most apparent feature of Figure 1 1 is the error bars showing the inter-quartile range of the 
length of the season. The contrast between NASS and the EDPM dates went to the edge of the 
anticipated differences due to disparities between compared datasets. We expected the 
variability in LoS to be driven by gradients in some climatic factors such as rain, temperatures, 
duration of daylight etc. What we found in NASS reports was that the states with more 
variability in seasonal precipitation (viz., Nebraska and the Dakotas) had more variability in 


576 phenological timing. The EDPM did not have the precipitation in the list of phenological 

577 controls [Kovalskyy and Henebry, 2011a] and, therefore, the vast difference between observed 

578 and estimated IQRs in LoS came as a result of limitations in number of factors considered as 

579 drivers of phenological timing. Moreover, the EDPM could not take into account the progress of 

580 agricultural work in spring as well as other anthropogenic factors affecting the development of 

581 crops. Nevertheless, the central tendencies were captured quite well for soybeans. The SoS 

582 delays in maize became the reason for underestimation of LoS for this crop. Yet, with all the 

583 shortcomings, the EDPM estimates of phenological dates for all crops and years managed to stay 

584 within the range of state reports from NASS. 

585 4. Discussion. 

586 Planned as a validation study this experiment took the form of a comparison between products 

587 while still providing insight on the performance of the EDPM +VegET scheme. In this context 

588 the discrepancies between estimates found in this study have to be considered just as relative 

589 indicators of better or worse perfonnance. Lacking the actual spatially explicit observations, we 

590 managed to obtain the reference points for the future application studies where the results will 

591 receive interpretation. It is clear now that the outcomes of this experiment helped reaching the 

592 goal of this investigation, yet they raised a number of other issues that need to be clarified. In 

593 each of the three sets of comparisons we presented spatial and temporal dynamics of error 

594 measures but we did not talk in details about the structure of uncertainties or about the reasons 

595 behind the observed patterns. Many of these issues are interconnected, and therefore we kept 

596 them for this section where the linkages can be explained. Every issue here is discussed in terms 

597 of its impact on the abilities of the EDPM alone and the EDPM plus VegET scheme to meet 

598 nominal performance expectations. We also present ideas about how these impacts can be 



599 mediated at this point and draw perspectives on possible corrections of the problems in future 

600 versions of the event driven phenology model. 

601 Comparison between the MODIS NDVI and the vegetation index produced by the EDPM had 

602 both temporal and spatial issues in performance. High r was definitely a plus to the EDPM, but 

603 the RMSE and bias of 2007 in prognostic mode pushed the perfonnance to the edge of what was 

604 expected of the model. Introduction of noise from the NDVI-TNDVI relationship could not be 

605 the reason for this error jump since such noise should have been present constantly and not just 

606 during late season drought on just about one-fifth of the study area. Apparently, the reaction of 

607 the EDPM to this development was too strong (residuals dropped to -0.25), most likely due to 

608 inability to account for irrigation. An appropriate solution for the 2007 error spike problem 

609 would be extra training of the EDPM on irrigated flux tower sites during the drought years. 

610 During other years, the bias appeared to be quite consistent throughout the area and could be 

611 arithmetically removed from the results. Possibly, the bias can be corrected by obtaining better 

612 estimates of background vegetation-free TNDVI values for growing season initiation as 

613 suggested by Zhang at a/., [2003]. This correction would, most probably, draw the overall 

614 RMSE close to 0.1 level. This performance mark was also achieved through the data 

615 assimilation. 

616 Patterns in temporal dynamics of residuals constitute a problem that cannot be corrected with a 

617 simple transfonnation. It requires collecting new data for parameterization of phenophase control 

618 module in the EDPM. Inclusion of precipitation as a control variable for phase transitions should 

619 help to address the issue of temporal variability in PTPs within states in addition to increasing 

620 the overall accuracy of the phenophase control procedures. With the current level of accuracy, 

621 we should refrain from interpreting the results based on uncorrected (prognostic) daily NDVI 



622 data in places where the variability of residuals goes beyond the level of two seasonal standard 

623 deviations. This warning, however, would be less applicable for time averaged (weekly or 

624 monthly) or composited prognoses. Meanwhile, the NDVI outcomes received from the EDPM's 

625 data assimilation scheme carried significantly smaller traces of phase control errors. Therefore 

626 further analysis can be conducted on the retrospective 1DKF corrected daily VI records and 

627 interpretations would be valid throughout the study region. The issues with temporal stability in 

628 perfonnance also came out in the ET a estimates produced by the EDPM+VegET scheme. 

629 Exceeding the expectation in r and RMSE in comparison with MOD 16 product, the results from 

630 prognostic mode exhibited a small bump and a dip of similar magnitude in temporal dynamics of 

631 the residuals. These fluctuations appeared exactly in the times of phenological transitions from 

632 green-up to reproductive phase and then from reproductive phase to senescence respectively. 

633 Present in results from all three years, the features indicated a systematic problem in 

634 phenological control module of the EDPM that, if removed, could further increase the 

635 perfonnance of the coupling scheme. In retrospective mode the results still had the issue of the 

636 early season overestimation. This indicates that while decreasing the level of variability in 

637 residuals the assimilation of MODIS NDVI could not completely suppress all the setbacks for 

638 the EDPM+VegET scheme. Improvements in the functioning and parameterization of 

639 phenological phase control module requires further training on long tenn flux tower records that 

640 will be undertaken in the future. However, all observed magnitudes of the deviations in temporal 

641 pattern would not pose a significant obstacle for the use of these results in further analyses. 

642 From the comparison of the EDPM+VegET scheme outcomes with ET a estimates from Mosaic 

643 LSM, we received the diverse spatial dynamics in r and RMSE complemented with clear 

644 seasonal patterns in temporal dynamics of residuals. These discrepancies persisted even after the 



645 assimilation of MODIS data into the EDPM and VegET results. In fact, the pattern became even 

646 more pronounced since the variability in residuals dropped. It is most likely that a better 

647 sensitivity of the EDPM to ongoing weather conditions contributed to the temporal dynamics of 

648 differences between two ET a estimates as the energy balance scheme in NASA’s Mosaic LSM 

649 [Koster and Suarez, 1996] uses static climatological trajectories of leaf area index as a 

650 phenology driven factor of canopy resistance. However, other patterns could not be explained 

651 entirely by the lack of sensitivity to contemporaneous vegetation development in the Mosaic 

652 model. It is also possible that numerous discrepancies came out as consequences of different 

653 assumptions about land cover on the 0.125-degree NLDAS grid and/or the ET flux partitioning 

654 between canopy and underlying soil. 

655 A unique feature of this study was the comparison of growing season metrics estimated by the 

656 EDPM with ones reported to NASS. In our analysis, we were missing proper geographic and 

657 temporal precision in the NASS reports for each of the five states. Nevertheless, we tried to 

658 preserve the temporal and spatial variability of growing season dates by organizing our SoS and 

659 EoS estimates to match the structure of reference data. We also kept in mind the fact that the 

660 transition points in NDVI dynamics and the actual phenological event for crops had different 

661 physical meanings. A good matching was achieved between reported and estimated state 

662 averaged SoS and EoS. Their variability, however, became problematic for the EDPM giving the 

663 ground to include more controlling variables into the automatic estimation of phenophase 

664 transition dates. 

665 Despite all the issues listed in this section the overall impression from the comparisons is quite 

666 positive for the VegET+EDPM coupling scheme. The scheme managed to keep the departures 

667 from references within nominal boundaries. The results matched and even exceeded most of the 



668 expected measures of model performance obtained on point based validations [Kovalskyy and 

669 Henebry, 2011a, 2011b]. The biggest problem for the TNDVI trajectories estimated by the 

670 EDPM was the model’s overreaction to late season drought in 2007 that accentuated the usually 

671 small underestimation. Meanwhile, the ET a estimates followed closely the reference records 

672 from MOD 16 products. Even in the worst cases, the error measures in ET a were also comparable 

673 with those of Senay [2008], Mu et al. [2007], and Abramowitz et al. [2008]. Remarkably, this 

674 level of perfonnance was achieved during the spatially explicit deployment of the coupled 

675 models. Plus, the results from the scheme were complemented with estimates of phenological 

676 metrics for grassland and crops that matched well the central tendencies of NASS reports. 

677 Combined with the ability of the scheme to produce daily estimates of vegetation index and 

678 actual evapotranspiration the perfonnance characteristics of the VegET+EDPM coupling scheme 

679 justified its use in a real life application study. 

680 The lessons learned from this experiment will help to analyze and interpret the results of the 

681 greater investigation of recent shifts in the phenology and ET regime in the Northern Great 

682 Plains. After the undertaken comparisons we can confidently say that consistency of received 

683 errors still allows for the trend analysis especially after correcting with MODIS observations. 

684 The delays of season starts in maize will be accounted for in the assessment of inter-annual 

685 variability of growing season parameters. Also, we intend to scale the variability in phenological 

686 dates from the EDPM to match the variability in NASS reports through inclusion of precipitation 

687 in the phenophase control mechanism. Special attention will be paid to the peripherals of the 

688 study region as those are most likely to carry land cover mapping errors. Finally, we will use 

689 appropriate testing methods and critical values when relating the shifts in ET a regime to crop 

690 cover change insuring a more conservative interpretation of their correlation. 



691 5. Conclusion 


692 The purpose of the experiment described in this paper was to provide the rationale for the use of 

693 the EDPM+VegET coupling scheme in a spatially explicit application. Such rationale was 

694 attained via assessing the performance of the scheme through comparison of modeled variables 

695 with reference data. First, we compared the image time series of vegetation index produced by 

696 the phenology model with MODIS NDVI derived from MCD43C4 product. The expectations of 

697 model perfonnance in producing seasonal NDVI trajectories were met yielding r" of 0.8 ±0.15 

698 and RMSE of 0.1 ±0.035 for the entire study area. Retrospective correction of canopy dynamics 

699 with MODIS NDVI brought the variability in errors closer to the 0.1 level. Estimation of 

700 growing season metrics by the EDPM matched the NASS reports with reasonable accuracy - up 

701 to 2 weeks of difference in key dates. The estimates of actual evapotranspiration produced by the 

702 coupled scheme were compared with ET a from NASA’s Mosaic model from NLDAS and with 

703 MOD 16 data from MODIS land product suite. In both comparisons, the expected r~=0.7 ±0.15 

704 and RMSE =1.4 ±0.5 mm per day were met by the coupling scheme working in retrospective 

705 mode using MODIS observations for correcting seasonal trajectories of canopy development. 

706 Minor issues of model perfonnance were encountered during this experiment as well. The 

707 EDPM produced trajectories of vegetation index biased towards underestimation but the bias was 

708 relatively uniform in space and time and therefore removable. Actual ET estimates from the 

709 VegET±EDPM were closer to MOD 16 product while producing greater differences with Mosaic 

710 LSM that also had persisting spatial and temporal patterns in them. While spatial patterns in 

711 differences could be attributed to distinct assumptions about land cover in Mosaic LSM [Mitchel 

712 et al., 2004], the seasonal profiles of differences between our estimates and reference data 

713 exhibited clear patterns driven by phenology. The impacts of these issues on performance of the 
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EDPM and the VegET models, however, were relatively small and therefore they could not pose 
an obstacle for the analysis and interpretation of the outcomes. In general, this study provided 
sufficient assurance that the interpretations of future results derived the the planned spatially 
explicit application study will be valid and sound, provided that the detected issues are properly 
addressed in the analysis. 
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Figure Al. Comparison of ET a from Mosaic LSM with the ET a produced by EDPM plus 
VegET coupling scheme deployed in (A) prognostic mode and (B) diagnostic mode 
involving 1DKF assimilation. (/) Coefficient of determination (r 2 ); (**) Root mean square 
error (mm per 8 days) . 
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Figure A2. Spatial distributions of residuals (A) ET a EDPM+vegET - ET a Mosaic (B)ET a edpm with 
iDKF+vegET- ET a Mosaic- (0 seasonal means of residuals (mm per 8 days) ; (**) standard 
deviations of residuals (mm per 8 days) . 
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